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ABSTRACT 

'x,/  Thia  paper  is  concerned  with  a  singular  perturbation  analysis  of  the  two- 
dimensional  steady  state  semiconductor  equations  and  of  the  usual  finite 
difference  scheae  consisting  of  the  five  point  discretisation  of  Poisson's 
equation  and  of  the  Scharfetter-Gumoel  discretisation  of  the  continuity 
equations#  By  appropriate  scaling7*^ transform  the  sesd  conductor  equations 
into  a  singularly  perturbed  elliptic  system  with  nonsaooth  data#  The  singular 
perturbation  paraaster  is  defined  as  the  ainiaal  Debeye- length  of  the  device 
under  consideration.  Singular  perturbation  theory  allows  to  distinguish 
between  the  regions  of  strong  and  of  weak  variation  of  solutions,  so  called 
layers  and  saooth  regions,  and  to  describe  solutions  qualitatively  in  these 
regions.  This  information  is  used  to  analyse  the  stability  and  convergence  of 
the  discretisation  scheae.  Particular  emphasis  is  put  on  the  construction  of 
efficient  grids.  It  is  shown  that  the  Scharfetter-GuaMl  asthod  is  uniformly 
convergent,  i.e.  the  global  error  contribution  coming  from  the  continuity 
equations  is  saall  when  the  aaxiaal  aash  size  is  saall,  independent  of  the 
gradient  of  the  solution.  Layer  jumps  are  automatically  resolved.  The  five 
point  scheae  however  is  not  uniformly  convergent. j  Large  gradients  of 
f soldtlonsrequire  a  graded  aesh  if  solutions  Inside  the  layers  are  to  be 
(^resolved  accurately.  T3»is_  c§n  J.ead  to  an  Intolerably  large  number  of 
gridpoints.^  Therefore,  ji^pc9*nt  a  modification  of  the  five  point  scheae 
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SIGNIFICANCE  AND  EXPLANATION 


It  is  wall  known  that  potential  and  carrier  diatribationa  in  a 
a eel conductor  device  can  be  mathematically  deacrlbed  aa  a  system  of  three 
elliptic  partial  differential  equationa  aubject  to  mixed  Meounn  and  Olrlehlet 
boundary  condi tiona.  The  dependent  variablea  are  the  electroatatlc  potential 
and  the  hole  and  carrier  densities.  We  ahow,  after  scaling  the  variablea 
appropriately/  that  thia  ayaten  ia  singularly  perturbed)  thia  meana  that 
second  order  partial  derivatives  are  multiplied  by  a  small  parameter.  The 
physical  interpretation  of  thia  parameter  is  the  Debeye-length  ^  of  the 
semi conductor  device  under  consideration. 

Asymptotic  analysis  (for  ♦  0)  allows  to  make  qualitative  statements 
about  the  solutions.  Xn  particular  regions  of  fast  variation  ('layers')  and 
regions  of  slow  variation  of  the  solutions  can  be  identified  and  solutions  can 
be  separately  described  in  both  types  of  regions. 

We  use  this  information  to  investigate  the  properties  of  a  widely  used 
discretisation  scheme  of  the  semiconductor  equations.  Particular  emphasis  is 
given  to  the  construction  of  grids  which  allow  an  efficient  numerical  solution 
of  the  discretized  problem. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  KRC,  and  not  with  the  authors  of  this  report. 
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A  SINGULAR  PERTURBATION  APPROACH  FOR  THE  ANALYSIS 
OP  THE  FUNDAMENTAL  SEMICONDUCTOR  EQUATIONS 


Peter  A.  Markovich*,  Christian  A.  Ringhofsr  and  Siegfried  Selberherr* 

1.  INTRODUCTION 

In  this  paper  we  present  a  singular  perturbation  analysis  of  the  two- 
dimensional  steady  state  semiconductor  equations  and  of  the  finite  difference 
method  used  by  Frans  et  al.  (1982).  The  singular  perturbation  approach  works 
as  follows.  The  carrier  densities,  the  doping  profile  and  the  independent 
variables  are  scaled  to  (maximally)  0(1)  such  that  Poisson's  equation 
assumes  the  form 

(1.1)  X2A*  -  P  ,  (x,y)  e  a 

s  s  s 

where  ♦  ,  p  is  the  scaled  potential  and  space  charge  resp.  and  R  is  a 

8  8  8 

2 

domain  in  R  of  diameter  0(1)  representing  the  device  geometry  after 
scaling.  X  is  the  dimensionless  minimal  Debeye  length  which  is  small  if  the 
maximum  of  the  absolute  value  of  the  doping  profile  is  large.  This  is  the 
usual  situation  for  modern  devices.  Therefore,  equation  (1.1)  subjected  to 
mixed  Dirichlet-Neumann  boundary  conditions  and  supplemented  by  the  scaled 
continuity  equations,  represents  a  singularly  perturbed  elliptic  boundary 
value  problem  (see  Fife  (1973))  which  can  be  analysed  by  adapting  well-known 
asymptotic  methods  (like  matched  asymptotic  expansions).  This  was  done  for 
the  one-dimensional  static  semiconductor  equations  by  D.  Smith  (1980), 
Markovich  et  al.  (1982a, b),  Vasilieva  and  Stelmakh  (1977).  It  turns  out  that 
in  every  closed  subset  of  R#,  where  the  doping  profile  varies  'moderately* , 
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the  solutions  of  the  seal conductor  equations  are  approximated  uniformly  up  to 
0(X)  by  smooth,  slowly  varying  functions  which  are  independent  of  X  and 
which  fulfill  the  'reduced'  equations  obtained  by  setting  X  “  0  in  (1*1). 
Between  these  subsets  there  is  a  curve  across  which  the  reduced  solutions  have 
a  jump- discontinuity.  Physically,  this  curve  is  the  junction  between 
differently  doped  regions  of  the  device.  We  derive  equations  for  the  Unite 
of  the  reduced  solutions  as  the  independent  variables  tend  to  the  junction 
from  both  sides  and  show  that  the  jumps  of  the  reduced  carrier  densities 
depend  exponentially  on  the  potential  drop  across  the  junction. 

Close  to  this  junction,  that  is  in  sets  where  the  doping  profile  varies 
strongly,  there  are  thin  regions  (of  width  0 ( X| £nX| ) )  of  rapid  variation  of 
the  potential  and  the  carrier  densities,  so  called  internal  layers,  within 
these  layers,  the  solutions  are  qualitatively  and  quantitatively  described  by 
the  solution  of  the  layer  equation,  which  is  a  second  order  ordinary 
differential  equation.  The  i-th  derivatives  in  perpendicular  direction  to  the 
junction  are  of  the  order  of  magnitude  X 

The  analysis  Bhows  that  even  large  changes  of  the  doping  profile  within 
layer  regions  only  cause  0(X)  changes  of  the  solution  outside  the  layer 
regions.  This  property  of  the  semiconductor  equations  carries  over  to  the 
discretisation  scheme.  It  causes  discretisation  errors  occurring  in  the 
layers  to  decay  rapidly. 

We  also  show  that  the  electron  and  hole  current  density  component,  which 
is  perpendicular  to  the  junction,  does  not  exhibit  layer  behaviour,  while 
layers  may  very  well  occur  in  the  tangential  components. 

Boundary  layers  can  occur  where  the  reduced  solutions  do  not  fulfill  the 
boundary  conditions  asymptotically  (as  X  0+).  This  happens  for  example  at 


oxide- semiconductor  interfaces  (for  MOS-transistors )  end  as  Schodsky  contacts 
but  not  at  Ohmic  contacts  and  Isolating  boundaries* 

Using  the  qualitative  and  quantitative  information  on  the  solutions  of 
the  semiconductor  equations  we  analyse  the  widely  used  difference  scheme  which 
is  obtained  by  discretising  the  Laplace  operator  by  the  usual  five  point 
formula  and  by  applying  the  Scharfetter-Guamel  (1969)  discretisation  to  the 
continuity  equations.  Due  to  the  strongly  different  behavior  of  the  potential 
and  the  carrier  densities  inside  and  outside  layer  regions  it  is  apparent  that 
the  construction  of  grids  has  to  be  done  with  particular  care. 

He  demonstrate/  that  the  chosen  discretisation  of  the  continuity 
equations  is  uniformly  convergent/  which  means  that  for  every  grid  6  the 
global  discretisation  error  e(G)  fulfills  the  estimate 
(1.2)  e(G)  <  const. (h  +  k  +  X|  in  X| ) 

where  h,k  are  the  maximal  mesh  sixes  in  x  and  y  directions  resp.  and  the 
constant  in  (1.2)  is  independent  of  the  grid  and  of  X.  The  Scharfetter- 
Gummel  scheme  resolves  layer- jumps  accurately,  even  without  using  a  fine  grid 
inside  layer  regions. 

Contrary  to  this,  the  five-point  discretisation  of  Poisson's  equation  is 
not  uniformly  convergent.  The  linearised  scheme  is  uniformly  stable  (i.e.  has 
an  inverse  which  is  bounded  independently  of  the  grid  and  of  X),  but  large 
discretisation  errors  within  layer  regions  can  destroy  uniform  convergence, 
particularly  when  0(X)-mesh  sixes  are  chosen. 

Therefore,  in  order  to  achieve  a  certain  given  (global)  error  tolerance, 
it  is  necessary  and  sufficient  to  control  the  grid  (only)  for  Poisson's 
equation,  since  the  error-contribution  from  the  continuity  equations  only 
depends  on  the  maximal  grid  sixes  and  on  X. 


We  show  that  there  are  two  possibilities  of  grid-control  for  the  five- 
point-formula.  The  first  Is  s  lsysr-ignoring  grid.  That  asans,  all  grid 
sisas  ars  choasn  to  ba  ouch  larger  than  X  which  implies  that  only  vary  few 
Msh  points  ars  located  within  layer  regions.  We  show,  that  for  such  a  mesh, 
the  solutions  of  the  discretisation  scheM  of  the  semiconductor  equations 
converge  to  the  reduced  solution  and  we  give  an  error  estimate  for  this  case. 

Of  course  the  choice  of  such  a  mesh  only  makes  sense  if  one  is  not 
interested  in  the  solutions  within  layer  regions. 

Therefore  we  also  derive  a  layer-resolving  mesh,  obtained  by 
equidistributing  the  local  discretization  error  (see  Markowich  and  Ringhofer 
(1982c),  Ascher  and  Weiss  (1981,  1982)).  The  construction  of  this  Msh  is 
based  upon  the  fact,  that  the  global  error  of  the  scheme  is  less  or  equal  than 
the  (linear)  stability  constant  times  the  maximal  local  discretisation  error. 
This  requires  the  information  on  the  exact  solution  acquired  by  the  singular 
perturbation  analysis. 

The  so  obtained  grid  is  coarse  in  regions  where  the  solution  is 
approximated  by  the  reduced  solution  and  it  is  fine  within  layers  in  order  to 
balance  large  derivatives  of  solutions. 

Storage  restrictions  normally  allow  to  use  the  equidistributing  Msh  if 
only  vertical  or  horizontal  junctions  occur,  however  junctions,  which  are  not 
parallel  to  the  x-  or  y-axis  usually  require  too  many  grid  points.  In  the 
latter  case  a  rigorous  equidistribution  is  virtually  impossible.  The  reason 
for  this  is,  that  in  the  case  of  a  horizontal  or  vertical  junction  only  the 
y  or  reap,  x  derivative  is  large  (the  perpendicular  derivative  is  large, 
not  the  tangential),  and  therefore  only  the  y  or  resp.  x-grid  sizes  have  to  be 
chosen  smII  compared  to  X  while  the  Msh  sizes  in  tangential  direction  may 
be  independent  of  X.  xf  the  junction  is  not  aligned  to  the  coordinate 


systmt,  all  partial  derivatives  can  gat  large,  than  tha  mesh  aisas  in  both 
directions  have  to  be  small  coatpared  to  X  within  tha  layer.  Nevertheless  it 
is  not  necessary  to  abandon  this  approach.  We  show,  that  even  if  too  few 
grid-points  are  placed  inside  such  a  layer  the  discrete  solutiona  of  the  five 
point  formula  are  qualitatively  correct  inside  layer  regions  and  they  are 
qualitatively  and  quantitatively  correct  in  'smooth'  regions. 

Since  this  situation  is  not  completely  satisfactory  we  present  a 
modification  of  the  five-point- formula,  which  is  uniformly  convergent  and 
therefore  does  not  require  any  grid-restrictions  (except  sufficiently  small 
h,k).  Grid  points  can  be  placed  wherever  the  solution  is  needed.  Again,  the 
asymptotic  results  on  the  solutions  are  heavily  used. 

We  also  give  an  analysis  of  the  "finite  boxes'  approach  (see  Franc  et  al. 
(1982),  which  allows  gridlines  to  terminate  outside  layer  regions.  "Missing' 
difference  quotients  are  approximated  by  interpolation.  We  show  that 
convergence  is  not  influenced. 

The  presented  analysis  demonstrates  the  power  of  the  singular 
perturbation  approach  in  obtaining  information  on  the  analytical  solutions  of 
the  semiconductor  equations  as  well  as  in  the  construction  and  analysis  of 
numerical  methods,  particularly  as  far  as  grid  construction  and  error 
estimates  are  concerned. 


2.  SINGULAR  PERTURBATION  AHALYglg 

Aa  shown  by  Van  Roosbroeck  (1950)  tha  equations  daacr thing  potential 
distribution,  carrier  and  current-distributions  in  a  a emi conductor  in  tha  two- 
dimensional  static  case  are 

(2.1)  edf  *  q(n  -  p  -  C)  (Poisson's  aquation) 

(2.2)  J  ■  -q(y  n  grad  ♦  -  D  grad  n)  (electron  currant  relation) 

n  n  n 

(2.3)  Jp  •  -q(y^p  grad  ♦  +  grad  p)  (hole  currant  relation) 

(2.4)  div  J  -  qR(n,p)  (electron-continuity  equation) 

n 

(2.5)  div  J  ■  -qR(n,p)  (hole-continuity  equation) 

P 

2 

for  (x,y)  e  (2  C  R  (where  Q  is  a  bounded,  convex  domain  representing  the 
device  geometry)  subject  to  Dirichlet  boundary  conditions  on  (Ohmic 

contacts)  and  homogeneous  Neumann  boundary  conditions  for  t,n,p  on  1*^ 
(isolating  boundaries)  with  30  *  T _  +  I\ .  Also  Dirichlet  boundary  conditions 

C  X 

for  J  =  J  +  J  and  oxide-semiconductor  interface  conditions  can  be  desired 
n  p 

for  the  simulation  of  certain  devices. 

However,  the  exact  formulation  of  the  boundary  conditions  is  not 
necessary  for  our  purposes  since  we  only  investigate  internal  layer 
phenomena.  The  numerical  and  analytical  treatment  of  boundary  layers  is 
completely  analgous. 

We  take  the  Shockley-Read- Hall  (SRH)-thermal  recombination  term 


(2.6) 


R(n,p) 


2 

np  -  ni 

t  (n  +  n  )  +  t  (p  +•  n. ) 
p  i  n  x 


In  order  to  model  high-injection  conditions,  the  SRH-term  has  to  be 
supplemented  by  more  complicated  generation-recombination  terms  (see  Schitt*  et 
al.  (1981)).  We  also  assume  the  validity  of  Einstein's  relation 


' hnifiMri 


For  simplicity  »•  wvm  that  °n'Dp  and  u  # w_  arc  constants.  Zn  'reality' 

•  ®  P 

they  are  weakly  varying  functions  of  n,p,  grad  f  and  of  the  doping  profile, 
nils  does  not  influence  the  following  singular  perturbation  analysis. 

An  existence  theorem  for  (2.1)-(2.5)  under  simplifying  assumptions  on  the 
boundary  data  and  on  the  device  geometry  is  given  by  Hock  (1974).  However,  no 
qualitative  information  on  the  solutions  can  be  obtained  from  this  theorem. 

Let  1  be  the  characteristic  length  of  the  device  under  investigation 
(for  example  the  length  of  a  diode). 

The  following  scaling  is  basic  for  the  singular  perturbation  approach 


(2.8) 


(2.9) 


♦  *  n  ■  "#  p  m 
a  s  r  •  - 


D  «  ~ 


IJ  U  ♦  ♦ 

_  n  „  p  .  n  .  p 

J„ - Jp  "  T*  *p  0 

s  D  qC  *s  D  qC  S  T  Fs  T 
n  p 


where  C  i»  max  |C(x,y)|  and 

(x,y)eQu8Q 


(2.10) 


x  y 

x.  mV  ys  t- 


Then  (2.1 )-(2.5)  reads#  after  dropping  the  subscript  s 


(2.11) 

(2.12) 

(2.13) 

(2.14) 


(2.15) 


X'  Lit  -  n  -  p  -  D(x,y,  X) 

J  -  -  n  grad  ♦  +  grad  n 


Jp  ■  -p  grad  -  grad  p 

div  “  5  S(n,p,YX) 
n  n 

div  J  *  -6  S(n,p,tX) 

P  P 


(2.16) 


_  X  2  eu„  „  iqn.  .2  .2 

k  -  (/)  --2^'  '.-dV-  "j-dV 

iqC  T  nnrpp 


(2.17) 


S(n,p,yX) 


.  »P  Z  7*X4 
n  +  p  +  2y2X2 


(2.11 )— (2.15)  holds  in  the  domain  A  “  {(x,y) | ( lx, £y)  e  Q}  and  is  subjactad 

to  (tha  scaled)  boundary  conditions.  is  the  minimal  Debaye  length.  Whan 

U  ,D  are  not  constant  D, D_  in  (2.9)  have  to  be  substituted  by 
n,p  n,p  n  p 

'characteristic'  values  5n,Dp  with  Dnp  -  0(5n  p).  The  scaled  continuity 
equations  (2.12),  (2.13)  have  to  be  changed  accordingly. 

In  the  sequel  we  also  need  the  scaled  Boltzmann  statistics 


2  2 

yV«  n. 


2  2 

-  Y2X2e  p 


Mow  assume  we  deal  with  a  silicon  device  with  characteristic  length 
l  ■  2.5  x  io  3 cm  and  the  doping  is  such  that  5  ■  1017cm-^.  Then,  at 
approximately  room  temperature  T  “  300K  we  compute  the  following  numerical 
values 

X2  -  0.4  x  io"6,  y2  ■  0.25,  B  •  B  "  0.25  . 

n  p 

Obviously  X  «  1  holds  while  the  other  constants  are  0(1).  Also,  the 
larger  the  maximal  doping  gets  (in  absolute  value),  the  smaller  X  gets.  For 
example,  C  -  10*vcm  J  gives  X  *  0.4  »  10  .  Therefore,  for  sufficiently 

large  doping,  we  can  regard  (2.11)-(2.15)  as  singularly  perturbed  elliptic 

Xk 

system  with  singular  perturbation  parameter  X  -  — . 

Markowich  et  al.  (1982b)  pointed  out  that  the  singular  perturbation 
analysis  requires  that  the  intrinsic  number  n^  is  much  less  than  the  maximal 
doping  5  and  that  C(x,y)/5  is  much  larger  than  X2  except  in  domains  of 
small  area  (i.e.  layers)  and,  of  course  that  X  is  sufficiently  small. 

In  order  to  analyze  internal  layer  phenomena  we  assume  that  the  scaled 
doping  profile  D  which  is  assumed  to  be  smooth  has  an  internal  layer  along 
the  y-axisi 


-8- 


j 


D(x,y,X)  •  D(x,y)  +  D(^,y)  +  0(X)  . 

A 

5,  D  have  a  jump  discontinuity  across  the  y-axis,  i.e. 

(2.18)  lim  D(x,y)  *  lim  D(x,y);  lim  D(?#y)  *  lim  o(rr,y) 

x+0+  x+0-  x-*0+  x+0- 

for  all  y  for  which  (0fy)  e  A  and 

*  .x  .  _C2^"X^ 

(2.19)  |D(f,y)l  <  cte 

holds  for  some  constants  >  0,  >  0  independent  of  x,y,X. 

Physically  this  represents  a  vertical  pn- junction  along  the  y-axis.  If 

A 

D  S  0  then  the  junction  is  abrupt,  if  (2.19)  holds  with  *  0  then  D  is 
assumed  to  be  continuous  and  the  junction  is  exponentially  graded. 

We  use  a  vertical  junction  because  this  heavily  simplifies  the  following 
analysis.  Later  on  we  state  the  generalization  to  curved  junctions. 

We  expect  the  internal  layer  in  the  doping  to  cause  layers  in  the 
dependent  variables  and  employ  the  'anaatz* 

(a)  *(x,y,X)  -  Wx,y)  +  *»(^,y)  +  ... 

(b)  n(x,y, X)  -  n(x,y)  +  n(^,y)  +  ... 

(2.20)  (c)  p(x,y , X)  *  p(x,y)  +  p(^.y)  +  ... 


(d)  Jn(x,y,X)  -  JR(x,y)  +  Jn(^#y)  +  ... 

(e)  J  (x,y,X)  -  J  (x,y)  +  J  (4,y)  +  ... 

P  P  Pv* 


where  the  dots  denote  a  power  series  in  X  starting  with  the  0(X)  term. 

The  functions  marked  with  '  -  ’  are  independent  of  X  and  denote  the 
reduced  solutions,  the  functions  marked  with  ' *•  are  the  layer  terms  which 
are  defined  for  t  -  ^  e  (-•,•)  and  all  y  in  the  device. 

The  layer  terms  are  supposed  to  decay  away  from  both  sides  of  the 
junction  j 

I 

1 
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(2.21 ) (a) 


♦(±«*#y)  ■  n(±«,y)  -  p(±“,y>  -  0 


A  A 

(2.21 )(b)  J  (±-,y)  -  J  (±“,y)  -  0  . 

n  p 

Now  we  insert  the  expansions  (2.20)  into  the  semiconductor  equations  (2.11)- 


(2.14). 


Neglecting  0 ( X)  terms  and  evaluating  away  from  the  junction  gives  the 
reduced  problem 


(2.22) 

(2.23) 

(2.24) 

(2.25) 

(2.26) 


0  *  n  -  p  -  D 

=•  -n  grad  +  grad  n 

J  =  -p  grad  -  grad  p  (x, y)  e  A,  x  *  0 
P 

div  J  *  B  S(n,p,0) 
n  n  r 

div  J  =  -0  S(n,p,0) 

P  P 


The  reduced  solutions  n,p,t|)  are  discontinuous  along  x  *  0  because  D  has 

*  . 

\  2 

a  jump  discontinuity  there.  In  the  sequel  we  denote  a  vector  a  ■  [  J  e  R  . 

ay 

Evaluating  close  to  the  junction,  but  to  the  left  of  x  =*  0  and  comparing 
0  terms  gives  the  left  layer  problem 


(a) 

♦  * 
tTT 

n  -  p  -  D 

(b) 

A 

n  ■ 

T 

A  A 

(n  +  n(0-,y))$T  -•  <  T  <  0 

(2.27) 

(c) 

A 

PT- 

A  A 

-(p  +  p(0-,y) ) 

(d) 

“x 

J  „  * 

0 

nT 

(e) 

a%. 

0 

PT 

where  we  set 

f(0±)  -  lira 

f  (x). 

Subscripts  denote  partial  derivatives. 

x+Oi 

Analogously  we  obtain  the  right  internal  layer  problem 
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(b)  nT  -  (n  +  n(0+#y))^T 

A  A  A 

(2.28)  (c)  pT  -  -(p  +  p(0+,y))*T  0  <  T  <  • 

(d)  JX„  -  0 


The  decay  of  the  layer  terms  and  (2.27),  (2.28) (d) , (•)  immediately  imply  that 

j(  *x 

(2.29)  OT  2  J  S  0 

n  p 

holds.  The  current  components  normal  to  the  junction  have  no  layers. 

Because  of  the  discontinuity  of  5  we  need  interface  conditions  at 

x  *  0.  These  conditions  are  obtained  from  the  (natural)  assumption  that  the 

solutions  ^,n,p,J  ,J  of  (2.1 1 )-(2. 15)  and  grad  ♦  are  continuous  along  the 
n  p 

y-axis.  This  implies  that  the  sum  of  the  reduced  and  layer  components  of  the 


expansions 

(2.20) 

have  to 

be 

continuous  < 

across  the 

y-axis 

(a) 

+(0-, 

-y) 

A 

+  ♦(0-,y)  “ 

♦(0+,y) 

+ 

*(0+,y) 

(b) 

n(0-, 

-Y> 

+  n(0-,y)  - 

n(0+,y) 

+ 

A 

n(0+,y) 

(c) 

p(0-, 

-y) 

A 

+  P(0-,y)  - 

p(0+,y) 

A 

p(0+,y) 

(2.30) 

(d) 

3f*<  o- 

n 

-#y) 

■  J*(0+,y) 
n 

(e) 

3y(o- 

n 

-*y> 

♦  Jy(0-,y) 
n 

-  5^0+ 

n 

,y) 

+  Jy(0+,y) 
n 

(f) 

^(O- 

p 

-,y> 

-  3£(0+,y) 

(9) 

^(0- 

p 

*#y) 

+  Jy(0-,y) 

-  0*10+ 
P 

,y ) 

+  Jy(0+,y) 
P 

and  by  comparing  o(j)  terms 


♦_<0-,y)  -  ♦_(0+,y) 


Integrating  (2.27),  (2.28 ) (b) , ( c)  gives 


(2.31) 


(2.32) 


n(T,y) 


f  _ 

n(0+,y) 

n(0-,y) 


(.♦<T'y) 

(eWT,y) 


D  , 
D  , 


p(t,y) 


f  p(0+,y) 


(e 


-♦<T,y> 


D  . 


I  p(0-,y)(e"^T,y)  -  1)  , 

and  (2.30) (a) ,(b) , (c)  give  the  interface  conditions  for 
densities 


T  >  0 

T  <  0 

T  >  0 

T  <  0 

the  reduced  carrier 


(2.33) 


U>  n<0+,y>  - 

(b)  p<o..y>  -  .Vo-.y>-«<o»,yi  j,  , 


The  layer  jumps  of  the  carrier  densities  depend  exponentially  on  the  voltage 
drop  accross  the  junction. 

The  reduced  problem  (2.22)-(2.26)  can  be  written  as  nonlinear  elliptic 
system  (after  a  lengthy,  but  simple  calculation) 


(2.34)  (a) 


Mi  +  — — — —  grad(2n  -  5)  *grad  ^ 

2n  -  5  2n  -  5 


8n  -  8  _  _ 

- - - -  S(n,n  -  5,0)  “  0 

2n  -  5 


(2.34)  (b)  An  -  (grad  n  -  ~ — —  grad(2n  -  5) )grad  ^i 

2n  -  5 


2n  -  D 


—  (AD  -  ( 8  -  0  )S(n,n  -  0,0))  •  0  S(n,n  -  D,0) 
S  n  p  n 


for  (x,y)  e  A,  x  4  o  subject  to  the  Neumann  and  Dirichlet  interface 


condition 


(2.34)  (c) 


5<0+,y)  -  •*0+'y,-*0-'y)  n(0-,y) 

(2.34)  (d)  n(0+,y)  -  5(0+,y)  -  a^<0"'y)_^<0+'y) (n(0-,y)  -  5(0-,y)) 


(2.34)  (•)  n(0+,y)*x<0+,y)  -  nx(0+#y)  -  n(0-,y) ^(O-^y)  -  nx(0-,y) 

(2.34)  (f)  (n(0+,y)  -  5<0+,y) ) *x<0+,y)  +  nx(0+,y)  -  5x(0+,y) 


-  (n(0-,y)  -  D(0-,y) M»x(0-,y)  +  nx(0-,y)  -  Dx(0-,y)  . 


(2.34) (•) , (f )  correspond  to  (2.3Q)(d) ,(f ) . 

The  reduced  problem  has  to  be  supplemented  by  boundary  conditions  (on 
3A)  derived  by  setting  X  ■  0  in  the  scaled  boundary  conditions  for  the 
full  problem  (2.11)- (2. 15)  on  Ohmic  contacts  and  isolating  boundaries  and  by 
making  boundary  layer  expansions  on  oxide-semi-conductor  interfaces  and 
Schofsky  contacts  similarly  to  (2.20). 

The  reduced  hole  density  p  can  be  computed  from  the  vanishing  space 
charge  conditions  (2.22). 

The  internal  layer  problem  is  obtained  from  (2.27),  (2.28)(a),  (2.31), 


(2.32) 

(2.35 

(2.35) 

(2.35) 

(2.35) 

(2.35) 

Markovich  et 


(a)  *TT  -  n(0+,y)e*<T'y)  -  p(0+,y)e"*{  T'y) 

A 

-  (5(0+, y)  +  D(T,y)),  T  >  0 

(b)  -  n(0-,y)e*(T'y)  -  p(0-,y)e"^  T,y) 

A 

-  (5(0-, y)  +  D( T,y) ) ,  t  <  0 

(c)  ♦(•/y)  •  ♦(-•vy)  -  o 

(d)  1>(0+,y)  -  1>(0-,y)  *  ♦(O-.y)  -  <<0+»y) 

(e)  iT(0+,y)  -  }T(0-#y)  . 

al  (1982b)  showed  that  the  ordinary  second  order  differential 
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equation  (2.35)  possesses  a  unique  solution  4  which  decays  exponentially  as 
t  ♦  *•  if  (2.18),  (2.19)  holds.  The  derivatives  fulfill 


(2.36) 


3^4  “Kjltl 

I — f  (T,y)|  <  K.e 
3T1 


where  ,K2  >  0  only  depend  on  i.  The  n  -and  p  layer  terns  have  to  be 
confuted  from  (2.31),  (2.32),  they  also  decay  exponentially  as  T  ♦  ±». 

After  having  shown  the  validity  of  the  expansion  (2.20)  (see  Markowich  et 
al.  (1982b)  for  the  one-dimensional  problem)  it  follows  that  the  solutions  of 
the  reduced  problem  (which  are  smooth,  slowly  varying  functions)  approximate 
the  solutions  of  the  semiconductor  equations  up  to  0 ( X)  outside  the  layer 
region,  which  is  a  strip  of  width  0 ( X | in  X| )  about  the  y-axis.  Within  this 
internal  layer  the  solutions  vary  rapidly,  for  example 


4<x,y,X)  - 


4(0+, y)  +  4(^»y)»  *  >  0 

4(0-, y)  +  4(4,y),  x  <  0 


There  the  solutions  are,  up  to  the  reduced  solution  evaluated  at  the  junction, 

x 

given  by  the  exponentially  decaying  (in  T  «  — )  layer  terms.  The  1-th 
derivative  in  perpendicular  direction  to  the  junction  fulfills 


(2.37) 


31  -i  -°2IXI 

I — 7  4<x,y,X)|  <  D  X  e 

3xx  1 


for  ,D2  >  0  independent  of  X  while  the  tangential  derivatives 
(y-direction)  are  0(1).  The  same  statement  holds  for  n  and  p,  but  not 

X  X 

for  the  'slow*  variables  J  ,  J  ,  which  do  not  exhibit  an  internal  layer. 

n  p 

They  are  -  even  close  to  the  y-axis  -  approximated  up  to  0(1)  by  the 

>x  — x 

corresponding  reduced  solutions 

So  far  our  analysis  did  not  give  any  information  on  the  tangential 

current  components  Jy ,  except  that  they  might  exhibit  layer  behavior. 

n  p 


A  complicated  asymptotic  analysis  (given  in  Markovich  (1982))/  which  usas 
tha  substitution 

A 

(2.38)  n  ■  aTu,  p  -  e  Tv 

shows  that 


(2.39) 


(2.40) 


>(T,y) 

^(T'y) 


.W0+,y)u^T,y)  _ 
.♦(0-,y)(#i(T,y)  _ 

V*°+'y>(e"*(T'y) 

a-W0-/y)(e-^T,y) 


u«y(o*y)» 

-  1)vy(0/y), 

-  1)vy(0,y), 


T  >  0 

T  <  0 

T  >  0 

T  <  0 


where  u,v  denote  the  reduced  solutions  u, v  (which  are  continuous  across 

x  *  0  because  of  (2.33)).  Therefore/  if  the  tangential  derivatives  of  the 

scaled  n  and  p  quasifermi  levels  +  ,  4  at  the  junction  do  not  vanish, 

n  p 

then  the  tangential  current  density  components  exhibit  layer  behavior. 

All  results  immediately  carry  over  to  a  horizontal  junction  by  exchanging 
the  x  and  y  coordinate,  however  a  curved  junction  extremely  complicates 
the  analysis.  Therefore  we  only  state  the  results,  the  proofs  can  be  found  in 


Markowich  (1982). 


Assume  that  there  is  a  curve  F  c  A  along  which  the  (scaled)  profile 
D  has  a  jump  discontinuity  (abrupt  junction)  or  decays  exponentially  (as  in 
(2.19).  Then,  for  points  (x,y)  sufficiently  close  to  the  junction  we  denote 
by  t(x,y)  the  closest  oriented  distance  to  T  (that  means  t  is  positive 
on  one  side  of  T  and  negative  on  the  other)  and  by  s(x,y)  ■ 

(s1 (x,y) ,s2(x,y) )  the  point  on  T  closest  to  (x,y)  (see  Figure  1). 
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Figure  1:  'Local  Coordinates' 


♦  ,tx(x,y)  \ 

For  (x,y)  e  T  we  call  n(x,y)  «  f  J  the  normal  vector  to  T  and 

-t  (x.v)  Vx'y) 

+  t  y  y  >  y  t(x  v) 

t(x,y)  ■  l  t  y ) '  the  tan9«ntial  vector  and  set  t  ■  •  It  turns  out 

that  the  asymptotic  expansions  (2.20)  remain  valid  when  x  is  substituted 

by  t(x,y)  and  y  by  s(x,y)  in  the  layer  terms.  That  means,  the  layer  is 

a  strip  of  width  0(A|An  A|)  about  T  (see  Figure  1).  The  exponential  decay 

of  the  layer  terms  occurs  in  perpendicular  direction  to  the  junction.  Both 

partial  derivatives  are  large  if  t  (x,y),  t  (x,y)  are  nonzero,  since: 

x  y 


(2.41) 


(a) 


(b) 


3x 


A  3t 


*  ^ 

3o  x 
1ft 


»tiyr.V  .  AmL  .  «  t  .1*! 111. 

3y  3y  A  3t  y  3s  x 


holds  where  ♦  and  fulfill  (2.36)  with  y  substituted  by  s.  The  first 

and  the  last  terms  on  the  right  hand  side  of  (2.41)  are  0(1).  Moreover 

-  ♦  •  ♦  > 

Jn*n,  Jp*n  *r*  continuous  across  the  junction,  while  Jn *t,  generally 

have  layers. 
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Th«  current  density  components  in  normal  direction  have  no  layer,  while 
the  tan9ential  components  may  very  well  have  one  (depending  on  whether  the 
tangential  derivative  of  the  quaaifermi  levels  vanish  or  not). 

Details  are  given  in  Markowich  (1982). 
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THE  BASIC  DISCRETIZATION  SCHEME 


For  the  discretisation  of  the  (scaled)  semiconductor  equations  (2.11)- 
(2.15)  ws  use  a  rectangular  grid  G  with  grid  sises  h^  in  x-direction  and 
kj  in  y- direct ion.  The  grid  points  have  coordinates  (x^.y^ )  where 

(3.1)  (.)  J).1  *  ')  *  *1 

holds.  The  subscripts  are  chosen  such  that  no  gridline  left  of  x  »  and 

below  y  ■  yfl  intersects  A  U  3A.  i  runs  from  0  to  N  and  j  from  0 
to  M. 

The  point  (x^,y^)  is  called  interior  gridpoint  if  its  four  neighbors 
(*i_1»yj>*  (xi+1*yj)'  <xi*yj+i>*  ^i'^j-l^  are  in  A  u  3A,  otherwise  it 

is  called  exterior  gridpoint. 

In  the  sequel  we  denote  the  evaluation  of  a  gridfunction  f  at 
(x^/y^)  by  fAj  and  we  define 

(3.1)  (b)  h  *  max  h.,  k  *  max  k  . 

i  j  3 

We  use  the  standard  five  point  formula  for  Poisson's  equation  (2.11) 


(3.2) 


h£j^Ol)]mn  -  p 

kj_1  JJ  "ij  Pij 


D(x^,y^,  A)  . 


For  the  discretisation  of  the  current  relations  and  the  continuity  equations 

we  need  the  function 

(3.2)  (b)  Y(z)  "  a  coth  s  . 

We  discretiEe  the  current  relations  (2.12)— (2.15)  by  the  following  difference 
scheme,  obtained  by  'exponential  interpolation': 


2 


2 


(3.2)  (c) 

(3.2)  (d) 

(3.2)  (e) 

(3.2)  (f) 

(3.2)  (g) 

(3.2)  (h) 


a*  -  nu  .  Vi.*  1  *i3  Vy  *  "lj 

nij  71  2  >  ht  ht  2 

jy  .  YA*in.r.  5.1)  1  nij .  jkJti-*  _V  n>.*»i  *  nu 

**  1  2  J  kj  ki  2 

jx  „  .Yf  jL±ti.“_V  .pu  _  *^,a.:Jy  fi+ia  *  pn 

pij  n  2  J  £  h4  2 


jy  .  .YfV,a*i.~  !ui  pAa+-u~  paj  .  "  jb  Pi>^  * pu 

pij  H  2  )  k  W  •> 


(jX  -  J*  ) 

'  n  n  ' 


h.  ♦  h  .  v’n  4  "n  )c.  +  k.  ...  .  ... 

1  *  »  *#j  1  >»j  j  j »  i»j  i»)”i 


(jy  -  J*  ) 

'  n,  ,  n,  j 


"  *nS,nij'pij'X) 


hi  +  hi-1  Pi,j  Pi-1,j  kj  +  kj-1  Pi,j  Pi,j-1 


(JX  -  ) 


(J7  -  Jy  ) 


"  ®pS<nij'Pij'X> 


(3.2)  only  holds  at  interior  gridpoints.  For  the  actual  computation  of 
solutions  (3.2)(c),(d)  and  (3.2)(e),(f)  are  inserted  into  (3.2)(g)  and 

(3.2)  (h)  reap. 

xt  is  an  easy  exercise  to  show  that  the  difference  schesw  (3.2) (c)-(h)  is 
equivalent  to  the  Scharfetter-Gusmel  discretisation  (1969),  however  the  form 
given  in  (3.2)  is  more  convenient  for  the  analysis.  (3.2)  is  used  by  many 
authors  (see,  for  example  selberherr  (1981)  and  Frans  et  al.  (1982)). 

The  boundary  conditions  are  discretised  in  the  normal  way  (see  Frans  et 
al.  (1982)). 
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In  the  following  Motions  we  apply  the  results  from  the  singular 
perturbation  analysis  to  construct  grids  which  allow  the  global  error  of  the 
discretisation  scheme  (3.2)  to  be  less  than  a  prescribed  tolerance  without 
employing  too  many  unnecessary  grid  points. 
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4.  AMhLYSIS  OF  TH>  SCHMOTTOUrgTOWtL  (SG)  METHOD 

Tha  on  of  tha  function  ylz)  as  a  factor  Multiplying  tha  diacrata 
derivatives  of  n  and  p  in  (3.2) (c)-(f )  appears  striking  at  tha  first 
glance.  Originally  tha  discretisations  of  tha  continuity  aquations  ware 
obtained  by  using  tha  exponential  dependence  of  the  corner  densities  on  the 
potential.  However#  by  applying  the  exponential  fitting  methods  described  in 
Section  5  for  Poisson's  equation  it  turns  out  that  (up  to  snail  terns)  the  80* 
method  is  the  only  discretisation  of  the  continuity  equation  whose  (discrete) 
solutions  have  the  aasui  asymptotic  properties  (as  1  ♦  0+)  as  the 
'continuous*  solutions  for  every  grid  with  R,R  suff  small  but  independent  of 
X. 

Taylor  series  expansion  shows  that 
(4.1)  Y(»)  -  1  +  0(s2)  as  s  ♦  0  . 

Therefore#  in  regions  where  differs  only  by  0(h  k)  from  its  four 

neighboring  values  *i#j-1  (for  sxample  outside 

layers),  the  discretisation  (3.2) (c)-(f )  is  up  to  0((h  +  k)2)  equal  to  the 
atandard-trapesoidal  rule  type-discretisation  of  the  current  relations  (with 
Y  =  1 )  However,  when  differs  significantly  from  one  its  four 

neighbors#  the  SG-method  behaves  differently.  To  demonstrate  this  we  apply 
the  method  to  the  one- dimensional  current  relation  and  continuity  equation 


(4.2) 

(a) 

n'  -  p'n  -  J  ,  -1  < 

n 

(4.2) 

(b) 

J'  »  0  ,  -1  <  X  <  1 

n 

(4.2) 

(c) 

n(-1)  -  n_  ,  n(1)  ■  n+ 

and  assume  that  the  potential  p  “  P(x,X)  is  a  given  function  with  an 
exponential  internal  layer  at  x  ■  0  of  width  0(X|fn  X|)  (see  Figure  2)  as 
given  by  the  singular  perturbation  analysis  of  section  2.  This  discouples  the 
semiconductor  equations. 
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Figure  2 t  Potential 


The  problem  (4.2)  can  be  thought  of  as  modelling  a  pn- junction  device  where 
the  junction  between  n  and  p  regions  is  at  x  -  0.  For  the  sake  of 
simplicity  the  recombination  rate  has  been  set  to  0  which  makes  sense  close 
to  thermal  equilibrium. 

A  simple  calculation  gives  the  exact  solution  of  (4.2) 


(4.3) 


+  J  e 
n 


1*x,X) 


X)ds 


-1 


(4.4) 


J 

n 


.*(1,X)  f  u 


-1 

At  first  we  investigate  the  SG-method  applied  to  the  homogeneous  initial  value 
problem 


(4.5) 

which  has  the  solution 

(4.6) 


n’  ■  ♦’n,  x 


n( x, X)  ■  n  e 


>  -1 ,  n(-1)  ■  n_ 

♦(x,X)-*(-1,A) 

e 
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The  SG- scheme  for  (4.5)  reads 


(4.7) 


-  *. 


i+1 


-  n. 


ri+1 


'  n 


i  i+1 


+  n. 


i  >  0»  n,  ■  n 


He  obtain  the  recursion 

(4.8) 


i+1 


with  the  growth  function 

(4.9)  a(s)  -  e2*  . 

Therefore  the  solution  n^  is  given  by 


(4.10) 


i-1 

n 

j-0 


V*0 

e  n 


♦(x^Xj-tf-l,!) 

e  n  . 

mm 


The  SG-method  integrates  the  initial  value  problem  (4.5)  exactly  on  every 
grid,  internal  layer  jumps  are  resolved  accurately. 

The  trapezoidal  rule,  that  is  (4.6)  with  Y  5  1,  gives 


(..in  Vi  ■  )"i 

with 

(4.12)  +(z>  -  —“7  . 

The  solution  of  the  trapezoidal  rule  is 


(4.13) 


i-1  "  % 

n  -  H  ♦(-3— 5 - i)n_  . 

j“0 


How  take  a  grid  which  is  such  that  no  gridpoint  is  placed  inside  the  layer. 

Let  Xj  be  the  largest  grid  point  'left  of'  the  layer,  then  is  already 

'right  of  the  layer. 
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Standard  theory  Implies  that 


♦(x  ,X)-*(-1,X) 

(4.14)  nj  "  n<Xj)  ■  e  n_ 

holds  if  h  is  sufficiently  small.  From  (4.13)  and  (4.14)  we  derive 


(4.15) 


1+1 


rI+1 


-  *. 


♦(x_,X)-*(-1,X)  * 

.  1  ♦ 


- 


lK  • 


♦(xI+i,X)-M-lfX) 

Since  n(xj+1 ,X)  «  e  n_  holds,  nj+j  onlY  approximates 

*1+1  "  *i  *txi+1»X)“WxI,X) 

n(xi+1,X)  if  $( - - - )  approximates  e  .  However 

$(z)  only  approximates  e2z  if  z  is  small.  Since  the  internal  layer  jump 

of  i>,  i.e.  <Mxi+j,X)  -  ♦(xI»X),  can  be  arbitrarily  large,  the  trapezoidal 

rule  generally  does  not  resolve  layer  jumps,  at  least  not  for  lay er- ignoring 

grids.  It  can  be  shown  that  the  trapezoidal  rule  converges  if  and  only  if  the 

grid  sizes  inside  the  layer  are  small  compared  to  X.  No  better  performance 

can  be  expected  for  the  boundary  value  problem  and  for  the  two  dimensional 

problem.  This  result  was  anticipated  by  J.  Barnes  (1976)  who  used  artyxsents 

based  on  physical  grounds.  The  superiority  of  the  SG-method  for  the  initial 

value  problem  is  due  to  the  exact  resolution  of  internal  layers  for  arbitrary 

grids. 

We  now  turn  to  the  boundary  value  problem  (4.4).  To  solve  the  recursion 

we  observe  that  J  3  A  holds.  Then  (4.4) (a)  has  to  be  solved  for  general 
"i 

A  by  using  (4.10)  and  variation  of  constants.  Then  A  and  the  second 
integration  constant  have  to  be  determined  from  the  boundary  conditions 


(4.4) (c).  We  obtain 


♦<X  .,X)-*(-1,X) 

(a)  \  m  *  n-  + 


(4.16) 


*  *  T'  “i 

+  I  Jn  h  • 

3-o  nj  3 


♦(xi#X)-*(xi+1,X)  t(x1+1,X)  -  ♦(x^,X)^ 

*(  9  ) 


where 


(4.17) 


„  _  ♦<1,X)-+<-1,X) 

n  ~  e  n. 


(4.16)  (b)  J 


n.  “  H-1  ♦(1,X)-Wx  ,X)  ♦(x  ,  X)-f(x  ,  X) 

1  I  .  J  '  h,.(  31  2 — 1 — 

j-0  3 


w(z)  *  e 


x  sinh  z 


holds.  The  structure  of  the  discrete  solution  is  analogous  to  the  structure 
of  the  continuous  solution  as  given  by  (4.3),  only  the  integrals  are 
approximated  by  sum. 

Standard  analysis  shows  that 


(4.18) 


|  /  '  e^^d.  -  "I  e"**Vw(2*lfi» 

Xj  j-I  3 


)|  <  Kh 


holds  if  the  layer  of  ♦  i»  not  in  [x ,x_).  Since  the  width  of  the  layer  is 

X  u 

0(X|tn  X|)  we  obtain  the  uniform  convergence  estimate 


(4.19) 


max  (|n  -  n(x  ,X)|  +  |J  -  J  | )  <  c(h  +  X| in  X| ) 

0<i<N  11  n  ni 


for  an  arbitrary  mesh.  If  the  mesh  sizes  outside  the  layer  are  constant,  then 
-  -2 

h  can  be  substituted  by  h  .  The  constant  c  only  depends  on  upper  bounds 
for  n, t  and  on  bounds  of  derivatives  of  n  and  ^  up  to  order  three 
outside  layers. 

The  estimate  (4.19)  carries  over  to  the  hole  continuity  equation  when 
n  is  substituted  by  p. 


Similar  results  hold  for  the  problem  with  nonzero  recombination  rate  and 
for  the  two  dimensional  problem.  We  do  not  present  the  proofs  here  since  they 
require  a  large  technical  apparatus. 

The  estimate  (4.19)  implies  that  the  error  contribution  of  the 
discretization  of  the  continuity  equations  only  depends  on  the  maximal  grid 
sizes  and  on  the  singular  perturbation  parameter,  but  not  on  the  mesh  inside 
the  layer. 

Actually,  this  estimate  can  be  refined,  such  that  the  right  hand  side  of 

(4.19)  is  substituted  by 

(4.20)  c(max  |t  |  +  Min  X| ) 

i.j  3 

where  Z ^  is  the  local  discretization  error  of  the  (two-dimensional)  SG- 
scheme  and  the  maximum  is  only  taken  outside  the  strip  of  width  0(X|tn  X| ) 
around  the  junction.  I only  depends  on  the  'local'  mesh  sizes 
h^,  h^_1 ,  k_j  ,k^_1  and  on  derivatives  up  to  order  four  of  n  and  (the 

local  discretization  error  is  computed  by  inserting  the  exact  solution  into 
the  SG  scheme  and  by  Taylor  expansion  (see  next  section)).  Therefore  the 
local  discretization  error  outside  the  layer  can  be  equidistributed  (see  next 
section)  and  the  grid  inside  the  layer  may  be  chosen  arbitrarily  at  least  as 
far  as  the  continuity  equations  are  concerned.  In  the  next  section  we  will 
show  that  the  grid  inside  the  layer  has  to  be  constructed  with  respect  to  the 
discretization  of  Poisson's  equation. 


5.  THE  FIVE- POINT  FORMULA 


The  scaled  Boltzmann  statistics 

(5.1)  n  -  Y2*2«*"*n»  P  -  y2*2#** ** 
transform  Poisson's  aquation  to 

(5.2)  X2A$  »  Y2X2(a  B  -  a  P  )  -  D(x,y,  X),  (x,y)  e  A  . 

Tha  linaarisation  of  (5.2)  with  raspact  to  ♦  along  soma  function  w  is 
givan  by 


2  >2 

(5.3)  X  Aw  -  Y  X  (a  ♦  a  p  )w  -  E(x.y.X),  (x,y)  «  A 

whara  E  has  a  layar  at  tha  sama  position  as  D.  The  function 


(5.4)  a  -  Y2X2(a  “  ♦  a  p  ) 

is  positive.  Sinca  discratisations  of  (5.3)  (including  boundary  data)  hava  to 
ba  solvad  in  every  stap  of  tha  (Mewton)  i tar at ion  of  tha  five  point  formula 
for  (5.2)  wa  lnvastigata  tha  five-point  formula  applied  to  the  linear 
singularly  perturbed  elliptic  problem 

(5.5)  (a)  X2Aw  “  a(x,y,X)w  -  E(x,y,X),  (x.y)  e  A 

subject  to  mixed  Neumann-Dirichlet  boundary  conditions  on  3A,  where  £ 
denotes  tha  outward  unit  normal  to  3A 


(5.5)  (b)  «|(M,  H|,M)  -0.  (>»)cU(J»)t-  >». 

c  1 

For  simplicity  wa  assume  that  tha  junction  is  abrupt,  i.a.  tha  function  B 

has  a  jump  discontinuity  along  soma  curve  rc  A  and  that  E  is  independent 

of  X.  Due  to  tha  built-in-potential  which  behaves  like  ±  to  y-  outside  a 

T  ^ 

pn-layer,  tha  function  a  is  uniformly  in  X  bounded  away  from  zero  outside 
tha  pn-layer.  For  simplicity  wa  assume  that  a  is  uniformly  bounded  away 
from  zero  everywhere,  that  means 
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t 


*  i 


(5.6)  0  <  o  <  a(x,y, X)  <  a,  (x,y)  e  A  U  3A 

with  a, a  independent  of  X.  This  assumption  does  not  influence  the  results, 
it  only  einplifiee  the  technical  approach  needed  for  the  analysis. 

Because  of  (5.4)  a(x,y,0)  has  a  discontinuity  along  T.  E  and 
<*(  •,*,<))  are  assumed  to  be  continuous  on  that  side  of  T  for  which  t  >  0 
holds  (one  sided  continuity).  The  five-point  formula  applied  to  (5.5) (a) 
reads 

(5.7) 

at  interior  grid  points  where  jW^^  is  the  discrete  Laplace  operator 


x  Lijwij  •  °ijwij  -  E(vV 


2 _ rritbil  wii .  lii  ivid) 

+  h,  .  1  h,  h,  .  > 


(5.8) 


LijWij  "  h 


i  "i-1 


i-1 


•  W,,  W,  ,  -  w. 


+ _ 2 _ _ lii  _  1LJ _ ZiiJcl) 

i  +  *j-i  Kj  Kj-i 


and  a 


ij 


“(Xj^.y^i). 


For  the  following  analysis  we  take  a  rectangular  device,  that  means  A 


is  the  interior  of  a  rectangle  in  the  x,y  plane.  Neumann  boundary 


conditions  can  be  discretized  in  the  obvious  way,  i.e. 


w_u~  >3  „ 


VW 


w  —  w 

—  tT - 4-1  *  VVyj>  and  analogously  for  uy.  This  leads  to  a  first 

order  approximation  of  the  boundary  conditions.  A  second  order  approximation 
can  be  obtained  by  ’mirror  imaging’,  that  is  by  introducing  the  exterior  grid 
points  (x_^ »y j ) »  (xp+i'yj)  for  the  approximation  of  ^(x^y^),  wx<xn'yj> 
and  analogously  for  wy.  After  arranging  the  gridpoints  row-wise  into  the 
grid  vector  w^^,  the  linear  system  representing  (5.7)  and  the  boundary 
conditions  can  be  written  as 
(5.9) 


Ito,!i<*,wh,k  "  fh,k 
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v-ATyTTT. 


is  equal 


where  ( f.  -  -  B(x. ,y.)  *t  interior  gridpoints  and  (f.  )  . 

ij  l  j  n#K  Ij 

to  the  boundary  data  at  boundary  gridpoints. 

Zn  Appendix  A  vs  provs  that  tha  system  of  difference  aquations  (5.9)  is 

uniformly  stable  in  1  in  tha  maximum  norm  for  every  grid,  that  means  there 

is  a  constant  L  independent  of  the  grid  and  of  X  such  that 

(5.10)  lti!k(X)l  <  L 

holds,  where  ve  denote  with  1*1  the  row  sum  norm  of  matrices  as  well  as  the 
maximum  norm  of  vectors. 

Stability  is  one  ingredient  for  convergence,  the  second  is  consistency. 
Therefore,  ve  insert  the  exact  solution  w,  that  is  the  solution  of  (5.5), 
into  the  difference  scheme.  Taylor  expansion  up  to  the  third  term  gives  the 
*  local  discretisation*  error 


“  X  VjjW(x^»yj )  -  «(xi,y^,X)w(xi,y^)  +  Etx^y^) 


(5.11) 


x2[“^7  *(5..^.) 


- r  wU,.  ,y.) 


'i-1 


3x3  3<hi  +  hi-1>  ax3  2i">  3(hi  +  hi-1> 


a  A  a  ki-i 

+  ^3  "(Vn1j)  3<X^  +  k^_1 )  "  ^3  *r<xi'T'2j)  3()Cj  +  k^) 

where  514/52i  e  *xi-i #xi) jni j ,n2j  ®  For  <5,11)  w®  assumed  that 

v  is  three  times  continuously  differentiable  in  A.  The  local  error 
component  from  the  first  order  discretisation  of  the  Neumann  conditions  at 
vertical  boundaries  is 


(5.’2)  *w-¥l*»'rtWV 

where  CQ  €  (xQ  ,x^ ) ,  C(J_1  «  ( 1 , x^ ) .  Similar  expressions  hold  at  horisontal 
boundaries.  For  general  grids  the  five  point  formula  is  of  first  order,  that 
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means 


(5.13)  |»lj|  -  0(h  +  k) 

holds  depending  on  bounds  of  the  second  and  third  derivatives  of  w.  For  a 
equidistant  grid,  i.e.  =  h,  =  k,  it  can  be  shown  by  taking  one  more 
term  of  the  Taylor  expansion  that  the  local  discretisation  error  is  of  second 
order  in  the  interior. 


He  denote  by  wex  the  vector  with  the  grid  values  of  the  exact 
solution  w  (row  wise  ordered)  and  get 


(5.14) 


Lh,k<X)(wex  “  Wh,k>  “  *h,k 


where  £.  .  is  the  vector  with  components  £,  .  The  stability  estimate 
n#k  ij 

(5.10)  gives  the  global  error  estimate 


(5.15) 


<Wex  -  whk*  <  L,*h,k'  * 


Therefore,  a  prescribed  global  error  tolerance  6  can  be  achieved  if  we 


require  that 
(5.16) 


'V  6 


holds  for  all  i,j  (L  is  given  explicitly  in  Appendix  A). 

Our  goal  is  now  to  determine  grid  sixes  ,  kj  such  that  the  local 
discretization  error  as  given  by  (5.11),  (5.12)  is  (approximately)  equal  to 
6.  This  process  is  called  equidistribution  of  the  local  error  (see  Markowich 
and  Ringhofer  (1982)).  The  obtained  equidistributing  grid  will  be  fine  inside 
layer  regions  (where  the  derivatives  vary  rapidly)  and  it  will  be  coarse 
outside  layers  (where  the  derivatives  are  0(1)). 

Equidistribution  is  normally  done  iteratively.  That  means  solutions  of 
the  discretization  scheme  are  computed  on  an  initial  mesh.  Then  the 
derivatives  in  the  local  error  are  approximated  for  this  initial  solution  and 
a  new  grid  is  determined  by  equidistribution.  A  second  approximation  to  the 
solution  is  computed  on  this  new  grid  and  so  on.  The  iteration  is  stopped  by 


an  appropriate  criterion,  e.g.  when  two  consecutive  solutions  differ 
insignificantly  or  when  a  solution  satisfies  the  local  error  condition  (5.16) 
on  the  aesh  it  was  computed.  For  the  nonlinear  problem  the  old  solution  is 
used  as  initial  guess  for  Newton's  method  for  the  computation  of  the  new 
solution  (continuation  method,  see  Ascher,  Christiansen  and  Russell  (1979)). 

We  will  now  show  that  the  singular  perturbation  analysis  can  be  used  to 
obtain  qualitative  and  quantitative  Information  on  the  equidistributing  grid 
inside  layers. 

Analogously  to  Section  2  we  derive  that  the  solution  of  (5.5)  has  the 

form 

(5.17)  w(x,y.  A)  -  5(x,y)  +  y(£lStgl ,s(x,y)  )  +  0(A) 

is  where  t,s  are  the  local  coordinates  introduced  in  Section  2.  The 
function  w  is  the  solution  of  the  reduced  problem  (5.5) 


(5.18) 


w(x,y) 


ZjXiXl 
a(x,y ,0 ) 


*  t 

which  is  discontinuous  along  r  and  w(-j,s)  is  the  interna 1- layer  term  which 
solves  the  layer  equation! 


(5.19) 

(a) 

-  «+(s)w,  T>0 

a  a 

(5.19) 

(b) 

wTT  -  o_(b)w,  t  <  0 

*  A  ^ 

(5.19) 

(c) 

w(0+,s)  -  w(0-,s)  -  w_(s)  -  w+(s) 

A  A 

(5.19) 

<d) 

*T<0+,s)  -  wT(0-,s) 

where  a+(s),a_(s) 

and 

w+(s),w_(s)  are  equal  to  o  and  w 

at  the  'right*  and  'left'  side  of  the  junction  T  resp.  and 


resp.  evaluated 
T  «  The 


layer  solution  w  then  is 


(5.19) 


(e)  w(t,s)  - 


\  x  >  0 
r-(«)«/*'U)  T,  t  <  o 


where 


(5.19) 


w  (s)  -  w^(s) 


It)  Y+(»)  - 


i  ♦/iSiil 

o-(e) 


(5.19) 


(g)  Y-(e)  - 


w+(«)  -  w_(e) 


1  ♦  /^iii 

o-(s) 


The  aeynptotic  expansion  (5.17)  holds  in  any  neighborhood  of  the  junction  F 
in  which  no  other  layer  occurs.  Me  remark  that  layers  can  only  be  due  to 
rapid  variation  of  E  and  to  boundary  conditions  which  are  not  fulfilled 
asymptotically  by  w  (boundary  layers).  Differentiation  as  in  (2.41)  gives 


(5.20)  (a)  w(x,y,X)  -  w(-  s(x,y))t  +  o(  -]■  -■■) 

3xA  X1  3T1  A  x  X1  1 

(5.20)  (b)  w(x,y , X)  -  ■—  s(x,y))t  +  °  ( — JZy) 

3y  X1  3tx  y  X 


inside  the  layer.  Therefore  the  local  discretization  error  (5.11)  fulfills 


h.  +  h  k.  f  +  k. 

(5.21)  |ttj|  <  const [ (— — f  ■-  |tx(xi»Yj) I  +  - — 1  |t  (x^y^)!  ) 


V  vs 

•exp(-  |t(x1,y^ ) | )  ♦  hi-t  +  hA  +  k +  k^] 


inside  the  layer,  that  is  for  |t(x,y)|  < —  | In  X| . 

/a 

Let  us  take  the  vertical  junction  x  «”o  at  first,  such  that  t  •  x, 
s  ■  ( 0 ,y ) .  (5.21)  now  reads 


r&. ,  ’  i. 


(5.22)  1*^1  <  const [■ 


h.  ♦  h 


✓5 


-  111  •xp(-  ~X  1**1 )  +  hi-1  +  \  *  J'j.,  + 


In  order  to  satisfy  the  equidistribution  condition  (S.16)  we  choose 

r* 


(5.23) 

(a) 

hi  “  C1  «p(-i  ixii ) 

(5.23) 

(b) 

"j  '  °2S 

for  |x^|  <—  X|in  X|  where  cf«c2  are  0  ( 1 )  constants.  Markowich  and 
/o 

Ringhofer  (1982)  showed  that  the  number  of  grid  points  in  x-direction  on  every 

gridline  is 

(5.24) 

inside  the  layer.  Due 

»U,.r  *  °<1> 

to  the  exponential  grading  of  h^  the  number  of 

gridpoints  is  independent  of  X.  The  number  of  y-gridlines  is 


(5.25) 


•ft)  * 


'layer  "'•6 
Therefore,  the  equidistributing  grid  requires 


(5.26) 


n  _  «*  Hy 

layer  layer  layer 


"(7> 


meshpoints  within  the  layer  in  order  to  achieve  an  0(6)  global  error.  The 
smallest  gridsice  hi  within  the  layer  is  0(X6)  and  the  largest  is  0(6). 
All  mesh  sizes  in  y>direction  can  be  chosen  as  0(6). 

Equidistribution  gives  a  comparable  number  of  gridpoints  for  horizontal 
junctions  and  for  horizontal  and  vertical  boundary  layers. 

We  now  assume  that  the  junction  T  is  not  aligned  to  the  coordinate 
system.  As  example  we  take  the  line 


(5.27) 


r  »  y  »  x 


such  that  t(x,y)  •  — j  (x  -  y)  holds. 


-33- 


(5.16)  and  (5.21)  give  bounds  for  the  mesh  sises 


(5.28) 


G  a 


(a)  ht  <  exp(— j  —  |xi  -  y ^  I ) 

1/2 

kj  <  c2X6  exp(-j  -j  \xL  -  y j  |  ) 


(b) 
2  X 


for  |x  -  y  |  <  ~  —  | In  X| .  The  number  of  meshpointe  inside  the  layer  is 
1  3  /2  /a 

at  least 

(5.29)  N,  -  0(-~r)  . 

iayer 


Since  the  x  and  y-derivatives  inside  the  layer  are  large  the  mesh  sises 

in  x  and  y  direction  depend  on  X  and  therefore  the  number  of  meahpoints 

inside  the  layer  increases  for  constant  6  as  X  decreases.  Even  assuming 

that  gridlines  terminate  outside  layer  regions,  this  leads  to  storage 

-3 

requirements  which  virtually  cannot  be  met.  Normally  X  <  10  ,  then  an 

accuracy  bound  6  -  10  implies  that  (in  the  order  of  magnitude)  10 
gridpoints  have  to  be  placed  inside  the  layer,  including  the  continuity 
equations  linear  systems  of  approximate  dimension  3  *  107  would  have  to  be 
solved  in  every  Newton  step. 

Therefore  equidistribution  has  to  be  relaxed  and  we  have  to  check  whether 
larger  mesh  sizes  also  lead  to  acceptable  discrete  solutions.  To  do  so  we 
eplit  the  matrix  L.  (X)  into 


(5.30) 


(1) 


Lh,k(X)  -  XK]l  * 


where  L^  ^  represents  the  discrete  Laplace  operator  (5.8)  (the  rows 

(2) 

corresponding  to  exterior  gridpoints  have  only  zero  entries)  and  L^  has 

diagonal  entries  at  rows  corresponding  to  interior  gridpoints,  it  has 

diagonal  entries  1  at  rows  corresponding  to  Dirichlet  boundary  points  and  the 

rows  corresponding  to  Neumann  boundary  points  have  two  nonvanishing  entries 

(2) 

representing  the  discrete  boundary  condition.  Obviously  L^  ^  is  nonsingular 


and 


(5.31) 

holds.  Also  we  obtain 

(5.32) 


*  <  max(1,  -^  +  h  +  k) 


where 


(5.33) 

holds. 

(5.34) 


h  *  min  h^. 


k  -  min  kj.  Assume  that 


Then  the  solution  of  (5.9)  fulfills 

(2)-1 


±\2 


h,k 


Lh,k  fh,k  +  °^h^ 


or  component  wise 

(5.35)  w±j  -  w(xi#yj)  +  0(h0  +  h„  +  kQ  +  k*)  +  o((£)2  +  (£)2) 

where  w  Is  the  reduced  solution  of  (5.5),  given  by  (5.18).  If  £  are 

n  k 

sufficiently  small,  the  discrete  solution  approximates  the  reduced  solution. 
Since  X  €  h^  X  «  k^  for  all  i,j  is  required  the  mesh  is  coarse  inside  the 
layer,  however  discretization  errors  which  occur  inside  layer  regions  do  not 
spread  out.  Of  course  layers  are  not  resolved  at  all  by  this  mesh.  By 
choosing  h  ■  0(6),  k  ■  0(6)  or  by  equidistributing  the  local  error  of  the 
continuity  equation  outside  layer  regions  with  error  bounds  the  reduced 
solution  of  the  semiconductor  equations  can  be  resolved  accurately  for  small 
X  on  such  a  grid.  The  bound  for  the  global  error  then  is 
0(6)  +  0 ( X)  ♦  0(|)2  +  (£)2). 

Therefore  even  a  constant  coarse  grid  h^  =  k^  £  h  >>  X  can  be  chosen  as 
initial  grid  for  the  iteration  procedure  and  the  computed  initial  solution 
will  be  close  to  the  reduced  solution. 
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He  now  investigate  the  behaviour  of  the  discrete  solutions  when  the 

condition  Z  snail  is  neglected, 
n  9c 

It  suffices  to  investigate  the  one- dimensional  constant  coefficient 
version  of  (5.5) 


(5.36) 


(a)  X  w"  -  au  +  E(x),  -1  <  x  <  1 


(5.36) 


(b)  w(— 1 )  -  w(-1),  w( 1 )  -  w(1) 


"15  . 

where  E  has  a  jump  discontinuity  at  x  *  0  and  a>0.  w  *  “  ~  “e 
reduced  solution.  Therefore  no  boundary  layers  occur,  only  an  internal  layer 
at  x  ■  0. 

The  one-dimensional  finite  difference  analogue  to  (5.7),  (5.8)  on  the 


constant  grid  with  h^  =  h  «  —  is 


(5.37) 


(5.37) 


(a)  ♦  <wi+1  -  2wi  +  wi_1 )  -  wt  -  E(Xi) 

(b)  wQ  *  w(-1),  wN  ■  w(1) 


(5.38) 


♦  -S 


He  regard  (5.37)  as  two-step  recursion.  The  characteristic  equation  is 


(5.39) 


$2(r2  -  2r  +  1)  -  or  -  0 


The  roots  of  (5.39)  are 


(5.40) 


(a)  r+( ♦) 


2 *2  +  q  +  />q»2  a2 


(5.40) 


(b)  r  (♦)  -  ^  °  ,  . 

2 


An  easy  calculation  gives 

(5.41)  (a)  lim  r  (♦)  -  +•,  lim  r+($)  -  1,  r+(#)  +  monotonically  as  ♦  ♦  • 

f*0+ 


(5.41)  (b)  lim  r  (♦)  -  0,  lira  r+($)  -  1,  r_($)  t  monotonically  as  ♦  ♦  w. 

♦>0+  "  fMi 


:  •  'wransBim 


The  difference  scheme  (5.37)  can  be  solved  explicitly 


i-H 


1  N-1  i-1 

(5.42)  w.  -  r*  *’a(N)  ♦  r‘b(N)  -  -  (  l  r*'3B(x  )  +  l  r  3B(x,)) 

/Zf77  3-1  3-°  '  1 


where 


(5.43) 


-  N~  <■»  «•  -S 

WN  ~  r-W0  W0  “  WNr+ 

•IMI  -  t  V  1’1111-  r  » 


-  -  (e)  '  -  (r> 


with 

(5.44) 


(a)  wQ  -  w(-1)  + 


1  r  -j 


/W  *  a2  3‘0 


L  r*  K<v 


(5.44) 


1  H  a 

(b)  w..  -  w(1)  + - - -  l  r_3S(x.) 

UA  "  J 


N 


/W  *  a2  *•» 


holds.  He  take  H  «  2L  such  that  x^  *  0  holds.  For  simplicity  we  assume 


(5.45) 


(a)  E(x)  - 


-1,  x  <  0 

1,  x  >  0 


Then,  using  (5.40),  (5.41)  we  obtain  the  asymptotics 


(5.45)  (b)  w±  -  <( 


i 

a 


(rt  -  1)/WT7 


r*"L  +  0(r"L)  +  O(r^),  i  <  h 


-I* 

a 


(1 


-  r  )/  4«$2  +  a2 


r*  **  +  0(r+L)  ♦  0(rS,  i  >  L 


The  first  term  on  the  right  hand  side  of  (5.45)  is  the  reduced  solution 


w(x)  • 


-,  x  <  0 


and  the  second  is  the  discrete  internal  layer  term 


-  x  >  0 


which  decays  away  from  x  ■  0  to  both  sides. 
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The  situation  for  the  two-dimensional  problem  is  completely  analogous. 
The  following  grid  construction  strategy  can  be  employed 

(A)  equidistribution  outside  layers 

(B)  'limited'  equidistribution  within  layers  ('limited'  means  that  not 
more  than  an  a  priori  defined  number  of  gridpoints  depending  on 
storage  restrictions  is  allowed  within  each  layer). 

A  discrete  solution  of  the  five  point  scheme,  which  qualitatively  agrees  with 
the  solution  of  the  continuous  problem  and  which  agrees  quantitatively  with 
the  continuous  solution  outside  layers  is  obtained  using  the  strategy  (A) , 
(B). 

This  however  requires  that  gridlines  are  allowd  to  terminate  outside 
layers.  We  now  discuss  the  situation  depicted  in  Figure  4. 


fj+l 


rj-l 


Xi-1  Xi  Xi+1  xi+2 


Figure  4t  A  terminating  line 

The  yj-gridline  terminates  at  x^,  the  point  (xi+1,yj)  is  no  grid  point. 

We  approximate  the  'missing*  x-dif ference  quotient  by  linear 
interpolation  between  the  (j  +  1)-st  and  ( j  -  1)-at  y  level. 
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(5.46)  ViiaVM  - _ Hiii.  + _ h _ wi+ij-r*iii-i 

k^  +  k^  1  hA 


The  right  hand  side  of  (5.46)  is  then  substituted  into  jw^ j  as  given 

(5.8).  The  local  error  introduced  by  (5.46)  fulfills 


h?  h2 


<r  * r  <‘i  ♦  Vi’  ^'‘5*  “i-i’) 


(5.47) 


•  «D?ul 


3  txi-l'Xi+11Xtyj-1'yj+1] 


where 

(5.48)  D*u  -  lu^l  +  lu^l  +  |uxyy|  . 

% 

Missing  /-derivatives  are  approxiaiated  analogously.  The  local  error 
contribution  then  is  obtained  by  interchanging  h^,  k^  in  (5.47)  and  d^u  is 
obtained  from  D*u  by  interchanging  x  and  y. 

(5.47)  implies  that  grid-size  ratio  restrictions  have  to  be  assumed  in 
order  to  get  consistency. 

For  a  terminating  ygridline  we  require 


(5.49) 


ii  <  N-) 

h  <  c#  v — 

hi  kj 


4  c 


and  for  a  terminating  x- gridline 


(5.50) 


hi 

r  * c' 

j 


*i-i 


4  C 


where  c  is  a  moderate  constant.  Then  (5.47)  simplifies  to 
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(5.51) 


ID*U» 


tXi-1 ,Xi+1 1 Xlyi-1 ,Yi+1 1 


and  analogously  for  a  terminating  x-gridline. 

The  local  error  introduced  by  a  terminating  line  is  of  first  order. 

Similarly  the  SG-scheme  for  the  continuity  equations  is  changed  at 

gridpoints  where  a  line  terminates. 

Generally,  lines  are  only  allowed  to  terminate  outside  layer  regions, 
x  y 

since  there  D^u,  D^u  have  moderate  values.  The  decision  if  and  where  a  x 
or  y-gridline  terminates  is  explained  in  Franc  et  al.  (1982). 

At  last  we  show  how  the  five  point  scheme  can  be  adjusted  such  that  the 
resulting  difference  scheme  is  uniformly  convergent.  We  modify  the  five  point 
scheme  as  follows: 


2ox 


~  w, 


•  W, 


20? 


— U-  rljUJ+ll.!*!  _  _ ILlIzI)} 

+  kJ  „  '  k.  It,  „ 


w, .  -  w. 


j  J-1 


3-1 


aiJwi3  -  E(xi,yj) 


where  o^,  will  be  determined  by  exponential  fitting  (see  Doolan,  Miller 

and  Schilders  (1980)). 

We  set 

(5.53)  h^  ■  pXA,  «  p^A 

and  since  our  previous  analysis  showed  that  only  mesh  sixes  which  are  of  the 
order  of  magnitude  of  A  can  destroy  convergence,  we  assume  that  p*, are 
independent  of  A  (later  on  we  will  get  rid  of  this  restriction). 

Zn  order  to  derive  a  necessary  criterion  for  o*  ,  to  make  the 

—  "  J 

scheme  (5.52)  uniformly  convergent  we  assume 
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(5.54) 


w^  -  wtx^.y^)  +  o(1)  as  X  ♦  0 
from  this  assumption,  inserting 
expansion  (5.17)  into  (5.52)  and  taking  the  limit  for  X  ♦  0  gives 


and  confute  ,  '*(  from  this  assumption.  Inserting  the  asymptotic 


20  Wl...  .,S..)-w(T  ,*  )  w(T  /S..)-w(T>  ,  A****) 

(5.55)  lim  [ - il-( - itU— - il-il - - k±L±-il-) 


,  .  X  ,  X 
x*°  »1  *  »11 


1-1 


,  2<1  ,*(ti.-i-h,*iii~'*(th,‘hi  .  ,(Tu-,u)--(ri.i-i-‘ii> 

ay  .  fly  1  ^  J 

rj  pj-i  *3“1 


-  lim  a. .w(t. .,s,  .) 


x+o 


ij  '  lj  ij 


t(x  ,y  ) 

where  ^  J  ,  s^  -  a*xi,yj  ^  hold*.  He  now  assume  that  the  five 

point  star  {  <xi-1  >  •  (x^y^ )  *<xi+1  'yj )  •  (xi'yj-i  >  •  <xi'yj+i * )  1b  •"tirely  on 

that  side  of  the  junction  for  which  t  >  0  holds. 


Since 


t(x.#y.) 


(5*56)  Vi.  j  "  T(Vi'yj) - IT*-  *  P?V  W'  ei  *  ‘Wi1 

and  similar  expressions  for  t  ,  t  ,  T.  .  .  hold,  we  set  using 

1“»,J  Irjrl 

(5. 19) (e) 

2o*  -^T)  'IW'l1  ,  ,  /“ij  °i-1txl,l,yil 

(5.57)  _  13 J'  (5 - - - ^  - - - ) 


X  ,  X 
0  +0 
'  i  pi-1 


i-1 


2«?.  -*T>  <5Wyj'  . .  «5-i  WV s 

*  — u —  f* - - - - - )  - 


f?  + 

pj  *>1 


«? 


<5i 


For 


0  <  c  <  1  we  set 
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, _ _ 

(5<s6)  (.)  i)  P^tx(«1.y1)i  -1  _  i  -  VjVyr^. 

'  X  ”  X  ' 


a *  - 


(P1  +  Pj-1)(1  "  c11)ali 


(5.58)  (b)  ij  exp(-/a^  pfcv(xi,y1>)  -  1  1  -  exp(/a^  *>![..,  *v<Xi'yj> ) 

«?  '  *1 


Consistency  requires  that  lia  o? .  “  lia  -  1  and  therefore  we 

*  -X  ..  3  J  v  3 
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obtain  by  Taylor  expansion 


(5.59) 


cl) m  t«,,i'yj)'  1  •  cij  ‘  *J(vV 


since  t2  +•  t2  «■  1  holds, 
x  y 

(5.58)  holds  for  all  i,j  for  which  the  five  point  star  centered  at 
(x^.y j)  is  on  that  side  of  the  junction  for  which  t  >  0  holds.  On  the 
other  side  /a^  has  to  be  substituted  by  -/a^.  in  order  to  understand  the 
influence  of  o*^,  o^  we  take  a  constant  grid  =  h,  k ^  =  k  such  that 
p*  5  p*  5  py  holds.  (5.58)  givss  for  t  <  0  and  t  >  0» 


(S.60) 


(a)  o*  ™o(u  )•  ■  o{v  ) 

1 '  ij  '  ij''  ^lj  '  ij' 


where  o(u)  ■ 


2'  uij 


a( xi #yj)  tx(xi*y.) )h 


^Q(xi,yi)  t  (xi#y,)k 


We  easily  get 


(5.61) 


J  9  *9  I  |i  I 

0(u)  "  1  ♦  0(u  ),  u  ♦  0,  c(u)  *  4u  e  1  '  as  |u|  ♦  • 


Therefore  the  scheae  (5.52)  behaves  (asymptotically)  like  the  aidpoint  rule  if 


h  <<  X,  k  <<  X.  Zf  h  >>  X,  k  >>  X  the  scheae  behaves  like 


*ij*ij  "  ^vV  +  °(#  3  ♦  •  3  ) 


(5.62) 
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where  e^  >  0,  f  ^  >0,  e^  ^  ♦  f  ^  >0,  The  reduced  solution  w  is  resolved 
(much  aore  accuretely  by  (5*60 )  than  by  the  standard  five  point  scheme. 

These  results  hold  locally  for  the  variable  mesh  sise  case  (5,58). 

The  modified  scheme  (5.52)  is  the  five-point- formula  analogue  to  the  8G- 
method. 

A  problem  which  still  has  to  be  resolved  is  the  modification  of  the  five 
point  scheme  at  stars  which  are  located  at  both  sides  of  the  junction  (see 
Figure  5). 


Junction 


Figure  5t  A  'Crossing*  Star 


In  fact,  it  is  for  the  abrupt  junction  formally  not  correct  to  write  down  the 
five  point  scheme  at  (x^,y^)  as  depicted  in  Figure  5  since  at  least  one 
second  order  partial  derivative  of  the  solution  is  not  continuous  at  the 
junction.  The  correct  approach  is  to  substitute  the  five  point  formula  at 
such  an  'interface*  point  by  a  discretized  version  of  the  interface  condition 


(5.63) 


grad(w) 


y  |t*0- 


at  (x1,y^). 

For  the  case  of  the  vertical  junction  x  ■  0  exponential  fitting  gives 


the  discrete  interface  condition 


(5.64) 


V  »  v 

L,j 


/a(0+,y ,)  /a(0-,y  ) 

1  -  exp( - j-J-  hj  1  -  ®xp( - 5—J-  h^,  ) 

assuming  that  x^  ■  0. 

Than  It  can  be  shown  that  the  modified  midpoint  rule  Is  uniformly 

convergent  and  that  the  global  error  is  0(h)  ♦  0(k)  +  0(A)  for  each  grid 

(see  Doolan,  Miller  and  Schilders  (1980)  for  one-dimensional  problem). 

The  disadvantage  of  the  modified  scheme  (5.52)  is  that  the  junction  T 

has  to  be  known  explicitly  since  o*  ,  of  depend  on  t„,t„.  For  most 

13  *  y 

practical  applications  however  the  doping  profile  is  graded  and  only  given  at 

t 

discrete  points.  Then,  since  )  is  parallel  to  the  direction  of  steepest 

y 

descent  (or  growth)  of  the  scaled  doping  profile  D(x,y,A),  the  approximation 


(5.65) 


f  *)  .  p 

't„'  Igrad  Dl 


can  be  used  for  the  (numerical)  determination  of  the  fitting  factors  (5.58). 


I 


I 


APPENDIX  A 


Let  bi j  denote  the  i-j-entry  of  the  matrix  fc(  A)  -  DA^  ^(A)  where 
D  is  the  diagonal  matrix  which  has  -1  at  every  entry  representing  an 
interior  gridpoint  and  1  at  every  entry  representing  a  boundary  point*  Then 


(a)  b  <  0,  m  *  n 

cm 


(b) 


1  ,  if  the  m-th  row  stands  for  a  Dirichlet 

boundary  point 

■  a  ,  if  the  m-th  row  stands  for  the  interior 
3  point  (xi#y^) 

0  ,  if  the  m-th  row  stands  for  a  Neumann 

boundary  point. 


(c)  B.  .  (A)  is  irreducible  . 

it  f  JC 


He  take  the  mesh  function 


a  + 


6  m 

?ij 


xi  + 


where  a,  0  are  chosen  such  that  .  >  c  >  0  and  max  I  .  I  -  1  where 

13  i,j  3 

i, j  runs  through  all  interior  and  boundary  points,  a,  6  only  depend  on  A. 


Let 

4.  .  denote  the  vector  with 

It  §  1C 

components 

(organised  row  wise).  Then 

“ij*u  • 

is  an  interior  point. 

{Bh,k(X)^h,k)ij 

♦t)  • 

if  (x^y^) 

is  a  Dirichlet  boundary  point 

i 

B  ' 

if  (x^y^) 

is  a  Neumann  boundary  point. 

From  Doolan,  Miller  and  Schilders 

(1980),  Appendix  A  we  derive  that 

(A1 ) 

•B*1.  (X)* 
n,K 

1 

1 

■i"  (Bh,k(X)  V^ij 

7*?  fa  *ij*ij  '*ij^ 

*#  J 

holds. 
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The  right  hand  aide  of  this  estimate  la  independent  of  X  and  of  the 
grid  and  can  be  taken  aa  the  atability  conatant  L  since  the  inverse  of  D 
has  norm  one. 

Assume  here  that  the  scaling  of  the  independent  variables  is  such  that 
the  rectangle  3A  has  the  corner  points  { ix,-£y) ,{-£X, iy)  and 

(ix,£y).  Then  a  simple  calculation  shows  that  the  constants  a,  0  can  be 
chosen  such  that 


L  -  2(1* 


♦  £y)  ♦ 


1 

■India) 


holds.  Mo  better  stability  constant  can  be  obtained  from  (A1 ) 
An  analogous  proof  also  holds  for  nonrectangular  domains. 
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